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Abstract — Standard model predictive control strategies imply 
the online computation of control inputs at each sampling 
instance, which traditionally limits this type of control scheme to 
systems with slow dynamics. This paper focuses on distributed 
model predictive control for large-scale systems comprised of 
interacting linear subsystems, where the online computations 
required for the control input can be distributed amongst 
them. A model predictive controller based on a distributed 
interior point method is derived, for which every subsystem 
in the network can compute stabilizing control inputs using 
distributed computations. We introduce local terminal sets and 
cost functions, which together satisfy distributed invariance 
conditions for the whole system, that guarantees stability of the 
closed-loop interconnected system. We show that the synthesis 
of both terminal sets and terminal cost functions can be done 
in a distributed framework. 

I. Introduction 

Model predictive control (MPC) is a well established 
method of process control that has proven to be useful in 
numerous industrial applications in the past decades. One of 
the advantages of MPC is that it can be applied to large scale 
systems, with a considerable number of states and inputs for 
which hard constraints are often required [18], [22]. 

MPC requires that the control input at each time step be 
calculated by the online solution of an optimization problem. 
As a result, one of the drawbacks of MPC as a control algo- 
rithm is the delay introduced by the computation time that is 
imposed for the evaluation of functions, their first or second 
order derivatives, and for matrix operations, computations 
that are usually required for most optimization algorithms. 
This computational burden is also worsened when MPC 
is implemented for a large-scale plant of interconnected 
subsystems, case where the dimension of the MPC problem 
is multiplied by the number of subsystems. For certain 
industries for which the manufacturing process is slow in 
nature, this computational time is not an issue. 

However, multi-system applications have arisen where 
computing the input rapidly is essential for efficiency and 
stability. Control problems for networks of interconnected 
multi-agent systems such as traffic control [10], building anti- 
earthquake systems [17] , satellite formation flight [20], and 
wind turbine farms [19], have received plenty of interest in 
recent years. Due to the large number of inputs and outputs 
of this class of systems, distributed control is often required. 
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Efficient distributed optimization methods for solving such 
control problems can be found in [2], [11], [23], [24], 
[25]. From a practical viewpoint, such methods can be sped 
up by implementing stronger, more powerful computational 
hardware. Recent results in [3], [4], [7], [12], [23], [24], 
[27] however, have shown that by exploiting the special 
underlying structure of some MPC problems, the number of 
flops required for an algorithm can be reduced substantially, 
thus making MPC a more attractive solution for control 
problems where speed is essential. The authors in [4] propose 
an interior point method approach for solving the MPC 
problem in which they use a discrete-time Riccati recursion 
to solve the linear equations efficiently at each iterate. In 
[3] the authors propose a more efficient approach to linear 
algebra computations w.r.t the derivation given in [27]. 
Computational burden can be also overcome by distributing 
the necessary operations amongst different agents. To this 
purpose, the authors in [23] examine a distributed approach 
to optimal control problems and appropriate optimization 
methods. 

In this paper, we focus on extending these recent results 
on the computational time required for the control action for 
MPC problems with a special underlying structure arising in 
large-scale leader-follower systems, where the computational 
burden is distributed amongst the comprising subsystems, 
thus providing a certain independence that is usually required 
for these subsystems. In the first part of the paper, a stability 
analysis for leader-follower systems is presented, based on a 
linear feedback law, that allows us to construct local terminal 
sets and cost functions in a completely distributed way. 
Compared with the existing approaches based on an end 
point constraint, we reduce the conservatism by combining 
the underlying structure of the system with distributed opti- 
mization. This leads to a larger region of attraction for the 
controller. Then, we formulate a distributed MPC problem 
for this type of systems, using a terminal cost-terminal set 
approach and an efficient implementation of an interior-point 
algorithm using Mehrothra's predictor-corrector scheme for 
solving the corresponding optimization problem is presented. 
In particular, we show how the underlying Newton system 
can be solved in a distributed manner. 

The paper is organized as follows. In Section [EI] we 
present the formulation of the MPC problem corresponding 
to systems of the leader-follower type and then we investigate 
the stability issue for the current system in a distributed man- 
ner via a linear feedback law using a structured Lyapunov 
function approach. In Section ITl-BI we focus on decomposing 
the terminal state constraints required for stability as a 



Cartesian product using distributed set computations, after 
which we formulate the general centralized MPC problem. 
We then show how to restructure the original MPC problem 
in Section [Hi] as to provide computational benefits using a 
distributed version of an interior-point algorithm presented 
in Section |IV] 

II. Distributed MPC using the terminal-cost, 

TERMINAL SET APPROACH 

Large scale systems have attracted much interest from the 
control systems community in recent decades. In this paper, 
we focus on large scale systems of leader-follower type. The 
MPC problem associated with leader-follower systems can be 
found in a number of current applications such as platoons of 
vehicles [26], which is of great interest in the development of 
automated highway systems [10], or in the renewable energy 
industry such as the problem of controlling a wind turbine 
farm [19]. 

Platoon or leader-follower systems imply that each sub- 
system, from the second one onwards, is influenced by the 
previous. We consider linear time invariant systems, for 
which the dynamics of the first subsystem are: 

xl+^A^l+B 1 ^. (1) 

The dynamics for the remaining M — 1 subsystems are 
described by the following linear equations: 

x\ +1 = A l x\ + B l u\ + A^x]' 1 + B^u]- 1 , (2) 

where x\ G W Li and u\ G M. mi are the state and input vectors 
of subsystem i at time t, A 1 G R n * xn * and B i G ]R™* xm ' are 
the state and input dynamic matrices for subsystem i, while 

A i,i-1 g piXnn and B i,i-l g R n»xm»_ 1 are the matrices 

for the coupling dynamics which define the influence of 
subsystem i — 1 upon subsystem i. For these systems, we 
consider mixed state and input constraints of the following 
polyhedral form: 

Glxl + Glut < b\ (3) 

where G l x G R 9 * Xn \ G l u G R 9 > xm % the matrix [G* G l u ] G 
R 9iX ™ i+mi has full row rank and b l > 0. We employ stage 
cost functions for states and inputs of the quadratic formr|: 

l\x\,u\) = \^\x\f Qi + \\4f R ^ 

where Q 1 G W liXn ' and R} G W ruxmi are positive definite. 
For the stability analysis, we also express the dynamics for 
the entire system as follows: 

x t+ i = Ax t +Bu f , (4) 

where x t G R™ and u t G R m comprise the states and inputs 
of all the subsystems at time t and the matrices A and 
B are block banded matrices comprised of A 1 , A" -1 and 
B l , B l respectively. In a similar fashion we define the 
block diagonal matrices Qd and Rd comprised of Q l and 
K 1 , respectively. In order to ensure stability for the MPC 

'in this paper, we use the following notation: ||a:||p = x T Px. 



scheme that we define below, we use a terminal set-terminal 
cost approach [18], [22]. We define the following final stage 
cost of the form: 

4(x) = ||x|£ d , 

where matrix Pd G M™ xn is positive definite. In order to 
find Pd and also a terminal set Xf we search for a linear 
feedback law u t = KdX t , such that the system 

x f+1 = (A + BK d )x t (5) 

satisfies the following three properties [18]: 

A.l {(x,K d x) |x g X f } C {(x,u)|G> 4 + G[y < V } 

A.2 (A + BK d )x g X u Vx g X f 

A3 if satisfies the following property: 

£ f ((A + BKd)x)-^(x) + x T KjRdKdX + x T Q d x < 0, 
Vx g Xf. (6) 

The centralized MPC scheme for the leader-follower system 
described by dynamics dTJ-d2j based on a terminal set- 
terminal cost approach, given an initial state x and prediction 
horizon N, is formulated as follows: 

M JV-1 

^v(x) =min V V e(x\,4)+tt(x N ) 

x,u * — ' * — » 

i=l t=0 

s.t: dynamics <[T} and (ffji (7) 

G l x x l t + G l u u\ < b\ x l =x\Vi = l,..., M, 
xjv g Xf. 

It is a well-known result [18] that the above MPC scheme, 
under assumptions A.1-A.3, stabilizes the system ®, with 
the optimal value of problem @> Vat(x), as a Lyapunov 
function. Keeping in line with the distributed nature of our 
system, the control law Kd, the final stage cost if and the 
terminal constraint set Xf need to be computed locally. In 
the following sections we develop a distributed synthesis 
procedure under such structural constraints. 

A. Terminal Cost 

For a locally computed Kd, we employ distributed control 
laws u l — K % x % for each subsystem, with K l g W"-**™ 1 
and the resulting control law for the entire system will 
then be u = KdX, where the matrix Kd = diag(i^) 
is block-diagonal. For the terminal stage cost, we define 

M 

i?f(x) = ^^tf(x % ), where terminal stage costs for each 

i=l 

subsystem are of the following quadratic form: 

ti(x i ) = ~\\x i \\ 2 pi ,Vi = 2,...,M, 

where the matrix P % G IR n i xn ; is positive definite, such that 
P d = diag(P J ). 

Due to the block-diagonal structure of matrices Pd, Qd 
and Rd, we can rewrite (|6) equivalently as the following 
inequality: 

M 

Vjix 1 ) + V}(x\x l - X ) < 0, Vx g Xf, 

i=2 



where the left hand side is a sum of local functions Vj that 
have the following form: 

V}(x 1 ) = (x 1 ) T ((A l ) T P 1 A l -P l +Q l + {K 1 ) T R 1 K 1 ') x 1 
Vj{x\x i - 1 )=[{x i ) T {x l ~ 1 ) T ] 1 P 



X" 

r.i-1 



,Vi>2, 
B i,i-i K i-i and 

\i\T pi Ai,i — 1 



where A i = A l + B l K\ A 1 '^ 1 = A 1 ' 1 " 1 
matrices P z are of the following form: 

~{A l ) T P l A l -P' l + Q' l + {K l ) T R l K l (A l ) T P l A 1 ' 

We can ensure inequality (O imposing the following dis- 
tributed structure (see also [11] for a similar approach): 

Vj(x 1 )<q 1 {x 1 ) (8) 
V}(x i ,x i - 1 ) <q\x\x l ~ l ), Vi = 2,...,M, xel f (9) 

such that: 



M 

q(x) = q 1 {x 1 ) + ^2q i (x i ,X i - 1 ) < 0, Vx G X t . 

i=2 



(10) 



We consider that the functions q 1 do not necessarily take 
negative values and have the following quadratic form: 

q 1 {x 1 ) = (x 1 ) T W 1 X 1 



q 1 (x l ,x l - 1 )=[(x i ) T {x l - l ) T ]W 
{W l )1 2 W} 



X" 

r.i-1 



where the matrices W l 



are symmet- 



12 "22 

ric. Clearly, <?(x) is also quadratic function and thus can 
be written as g(x) = x T Wx, for an appropriate matrix W 
defined below. We now define the following optimization 
problem: 



min t 



(ID 



s.t: ML l {P\K\W x ) 4 0, Vi = l,. 
PF rJ, 



where MF(-) refer to the matrix inequalities ((8} and ((9) and 
the matrix W has the following block tridiagonal structure: 



wi 2 



w 1 - 
{wh) T 






w? 2 

+ 1F 2 3 2 








^r 1 











It is straightforward to see that if (fl~T) has an optimal value 
t* < 0, ensuring that W < and subsequently (TTOb holds, 
then © is satisfied. Note that we do not require that matrices 
Wi to be negative semi-definite. On the contrary, positive or 
indefinite matrices allow local terminal costs to increase as 
long as the global cost still decrease. This approach reduces 
conservatism in deriving the matrices Pj, and K^, However, 
problem (fTTT i is in a form that cannot be solved efficiently 



since it is not a convex problem. Subsequently, we show 
that ( fTTT ) can be expressed as a sparse SDP that can be 
solved distributively. To this goal, we employ the following 
linearization [22]: 



P l = {S l )-\ K l = Y l G- 



(12) 



ir, 



12 



G + G T — f/l + W, 



22 



We also define the following matrices in order to make 
constraints of the optimization problem in the following 
theorem more compact notationally: 

G + G T ~S 1 + W 1 
\G + G T -S' + W\ x 

m 2 ) T 

~A 1 G + B 1 Y 1 ' 

\A l G + B i Y l A 
[Q^-G 

(R^Y 1 
G T 



G 1 

& = 

B 1 = 

B l = 



G 



o 



S 1 



S 1 0' 
0/0 
0/ 



s l = 




/ 
I 



n l I 

Note that the linearizations (fT2l have been employed 
under the assumption that all the subsystems have the same 
dimension for the states, i.e. rii = rij for all 

Theorem 1: If the following SDP 



s.t: 



(13) 
(14) 



mm t 

% ( ff]^0,Vi = l,...,M 
y-i.i-i = yi-i j \/i = 2,...,M 
W 4 rl, 

where W has the same structure as W, has a negative optimal 
value t* < 0, then © holds. 

Proof: From (TPfl i we observe that S l >- and /i J > 0, 
which in turn implies: 

(S i -G) T {S i )- 1 (S i -G) ^ 

{ l u i I-G) T - i I( t i i I-G))p0, 
M 

By adding W l to the previous inequalities, we get: 

'(S 1 )- 1 

' o XI 



G 1 4 



G T 


" 







G T _ 







G 


0" 







G 



W\ (15) 



For i — 2, • • • , M, using (Q3) and the equality constraints 
Y t ' i — Y 1 ^ 1 and by applying the Schur complement to 
(fl4l we obtain: 



(A i ) T P i A i - p i (i i ) T P i i i ' i - 1 

>' + {K l ) T R l K l + G- T G- 1 0' 




which is equivalent to (0 if we take: 



G- T 










G- 1 



W{ 2 G~ 
[{W{ 2 ) T wi 2 \ [ o 

To transform inequality (0 into a linear matrix inequality 
of type (O, we use the same linearizations and the proof 
follows similar steps as those previously presented. As a 
result, the SDP dT3l l is equivalent to problem (fTTT l. and for a 
negative optimal value r*, (0 is satisfied. ■ 
Note that the SDP problem (fT3ll-(fT4l> can be solved offline 
either using a sparse SDP solver or some distributed opti- 
mization algorithm [23]. Since we have imposed a diagonal 
structure on the controller Kd = diag(A' 1 ), it follows 
that the system matrix A + BKd has a block bidiagonal 
structure. If the optimal solution r* of the SDP is negative, 
then the matrix A + BKd is Schur (all the eigenvalues are 
strict inside the unit circle). It follows that all the matrices 
A 1 + B l K i are Schur. 

B. Terminal Set 

To complete the stability analysis for system 01, which 
implies properties A.l - A. 3, we need to complete the design 
procedure by the computation of a terminal set Xf C K n , 

M 

defined locally (as a Cartesian product) Xt = II XI and 
equipped with invariance properties. 

First let us define the set of admissible states associated 
to the constraints (0 an d the specific linear controller K t : 

X 1 = {x l : (Gl + G\K l )x l < b 1 }. 

leading via the Cartesian product to a set in W 1 : 

M 

x = nr. 

<=i 

Assumption 1: The origin is assumed to be an interior 
point of the set X. 

We introduce the following formal definition of positive 
invariance in view of its use in the practical construction of 
the terminal set Xf. 

Definiton 1: A set VL C X is called positive invariant for 
system (0 if x f G VI it holds that x t+ i e il for all t > 0. 
As a standard approach in the MPC design [18], the terminal 
set Xf C X needs to be positive invariant for the nominal 
linear time-invaraint dynamics (0. This is a standard prob- 
lem in set-theoretic control theory and there are a number of 
ways in which can be computed (see e.g [16], [6], [8]). 

Due to the distributed nature of our system, such a general 
terminal constraint set cannot be used due to the introduction 
of coupling constraints between the states of the subsystems. 
We need to explore the possibility of finding a terminal 
constraint set, which preserves the structure of a Cartesian 
product: 

M 

X t = UX}, (16) 

i—l 

This will further enable a distributed use of the terminal 
constraint sets X\ for each of the subsystems. Is worth 
mentioning that for general systems the construction of a 



terminal set in the form given above can be cumbersome 
in distributed settings (see e.g. [5] for such a construction). 
However, for a system x t+ i = Ax t , where A has a 
special block bidiagonal structure and the admissible set 

M 

is expressed as X = II X'\ the computation of such an 

i—l 

invariant set X t = IL^Xf can be simplified by exploiting 
these structural properties. 

Without loss of generality the matrix A will be considered 
to be of the following form: 



A = 



i.e block lower-bidiagonal. The developments in Subsection 
IH-Al point to the construction of a distributed linear controller 
which allow us to assume the stability of the unconstrained 
local closed-loop system x?+i = Ax t around the origin. By 
the block lower-bidiagonal structure it follows that the matrix 
A is Schur (i.e. 



'A 11 










A 21 


A 22 











A 32 


A 33 




















q fiM,M-l 


£MM 



A(A) < 1) and consequently through the 

block lower-bidiagonal form of A, all the matrices A 11 are 
also Schur, for all i = 1, • • • , M. 

The dynamics for the comprising subsystems are: 



b t+i 



A n xl 



! = A ll x\ + A 



Vi = 2, . . . , M. 



(17) 
(18) 



1) Construction of Xj: By taking into account that the 
first subsystem is stable and its dynamics are not perturbed 
by the other subsystems, the computation of Xf C X\ as 
a positive invariant set with respect to ([PTt can be done 
easily through standard methods for LTI nominal dynamics 
available in [6], [8]. 

We note also that the boundedness of the set X\ will 
ensure boundedness properties for the set Xj. 

Remark 1: If Xj C X\ is invariant with respect to (TTTb 
and £ int(A^) then aXj is invariant and € int(aiAj) 
for any scalar a > 0. More than that, if < a < 1, then 
aXj C X x . 

2) Completing the construction of Xj: For the subsys- 
tems i ~ 2, . . . , M we require a different treatment. If we 
denote A l - l ~ 1 x\~ 1 — w\, the dynamics for the remaining 
subsystems can be considered as: 



h t+i 



,M, 



(19) 



where can now be viewed as an unknown disturbance 
for this particular subsystem, where w\ is bounded, i.e 
w\ is in a set We denote by w(-) G Myy the 

sequence Wo, W\, . . . , Wk of disturbances from the admissible 
set M w = M')K € W, Vfc e N}. 

'in the case of the second subsystem i = 2, we have uij = A 21 x], and 
by taking into account that Xj is positive invariant, it can observed that 
w? is bounded, i.e w? £ W 2 , where W 2 = A 21 X}. 



Definiton 2: The set O C X is a robust positive invariant 
set for a system x t +\ — Axt + w t , if starting from 0, the 
evolution of the system remains in O for all w(-) G A4w- 
We observe that Xf can now be computed as a robust 
positive invariant set (RPI) for the subsystem with the index 
i > 2, by exploiting the contractiveness properties of A" 
and the existence of explicit bounds on w\. The practical 
construction of such RPI sets is standard in the literature, 
see for example the procedures in [6], [14], [15]. In the 
following such a constructive procedure will be denoted by 
XI = RPI(X\W l ). 

Proposition 1: Let X\ = RPI{X\W l ) be an invari- 
ant set with respect to (1181 . having the origin as interior 
point. There always exists a scalar < a < 1 such that 
aX\ = RPI(aX l ,aW l ) preserve the invariance properties 
and additionally aX\ C Xj. 

Proof: The proof is an immediate application of the 
Remark [T]and the scaling properties of the RPI sets detailed 
in [13]. ■ 

With these (robust) positive invariance and constraint 
satisfaction properties we are able to propose a constructive 
procedure for X\ in a iterative manner, starting from the first 
subsystem and leading to an invariant set in R™, as presented 
in the following algorithm: 

1) compute X] 

2) for i = 2 : M 

1. compute W l = A i < i ~ 1 X^ 1 ; 

2. compute X| = i?P/(X 4 ,W i ) ; 

3) find < a < 1 such that aX\ C X u Mi = 1, . . . , M 

Since for the leader-follower systems described in this 
paper the matrix A + BKj is block lower-bidiagonal as 
well, we can use the procedure described above to compute 

M 

a terminal set of the form Xr = II Xi that satisfies the 

i=l 1 

properties A.1-A.3. Note that the distributed MPC controller 
presented below results in a larger region of attraction 
compared to other MPC schemes based on an end point 
constraint [2]. An additional novelty of our approach consists 
in the fact that all the computations for the terminal set 
and cost can be carried out in a completely distributed way. 
Note that this strategy for constructing sets Xf can also be 
extended to the case where A is block lower triangular, i.e 
subsystem i > 2 is affected by subsystems — 1. 

In this case, the sets W" would be constructed as W l = 
W M © • • ■ © W*'* -1 , where by © we denote the Minkowski 
sum: A © B = {x + y\x G A, y G B} and W ,;j = A^Xj, 
j = 1. 

We can now reformulate the centralized MPC problem for 
the entire system (01 as following: 

M N-l 

^v(x) =mm V V ^ t) <4)+4(4r) (20) 

x,u A — » A — » 

1=1 t=0 

s.t: dynamics ((TJ and @, 

Glxl + G>j < b\ G l x l N <f,Vt = l,... 7 M 



where we assume that the terminal sets X{ constructed 
previously are polyhedra described by G l x l N < /', with 
/ 4 >0. ' 

III. Problem restructuring 

We now propose to reformulate problem (l20l as to obtain 
a more suitable structure. We define the intermediary stage 
variables for subsystem i as: 

xi = [(xj) T K) T ] T eM n s 

where = rii + rrii and t = 1, . . . , N — 1. Next, we define 
the general decision variable z G 1" for (f20b as follows: 

z = [(zT-.VT] T , 

M 

where n = A^n^ and 



i=l 



TiT 



z* = [K) T (x\) T ...(v^ N _ 1 ) T (^-i) T (^) 

Now, in accordance with the general decision variable 
as defined above and in order to create a more compact 
and ordered final structure, we need to define the following 
matrices: 



E l = [I ni 0] G W 1 ' 
A M-1 = B^- 1 
Q' = diag(Q 4 ,i? 1 ) G R n » xn » 
Q i = diag(ii i ,Q i ,...,Q i ,P i ), 



A 1 = [A 1 B 1 } G 

] G M™^" 1 - 1 



where Q l G R N ^ xNn - has N-l Q l blocks in its diagonal. 
Using the intermediary stage variable we can rewrite the 
equality constraints in < f20b for subsystem i as: 



E*xJ +1 = AX 



We now recast (|20] > as: 



min z T Hz 



s.t : Gz < b, Cz 



(21) 



(22) 



where H G M nxn is diag(Q 4 ), with i = l,...,M. We have 
included the equality constraints for each subsystem in (|20T i 
in Cz = c, where c G Ar ™ 1 and C G 

are the following: 



,4 2 



A 1 ^ 

0(Ar-l) ni 4 



0(iv-i 



- yl 2 a; 
)n 2 ,i 



Or 



'(N-l)n M A 



C = 















C 21 


C 22 













C 32 


G 33 





■ (23) 




















qm,m-i 


(jMM 





In (O the matrices C" £ I 
have the following structure: 



>NmxNri 



C" = 



-B l 







E* 
-A* 






E 4 

-A* 





E' 



for i = 1, . 







,M, 



-A 1 



whilst the matrices C 1 ' 1 " 1 G R^x^-i, fori 




= 2. 



have the following structure: 

- _ B i,i-i 


















A i,i-1 









_ A M- 



The inequality constraints in (T20l > have been recast as Gz < 



M 

b, where 6 £ M q , with q = ^(Nqi + q) and G € 

i=l 

have the following structure: 

&=[(bT,-.-,(b AATlT 

where 



]>qxn 



N — l times 



b'=[(i'-GX) T ,(t'f 



and 



G 



G 1 









G 2 



m T ,(n T ] T , 



whose matrix blocks G l £ 
\Gi . 





mNqi+qxNlli 







(24) 



are: 



G' 











u 



















G i 



IV. Primal-dual interior point method 

Primal-dual interior point methods are very efficient op- 
timization methods which employ the Karush-Kuhn-Tucker 
(KKT) conditions, that are both necessary and sufficient 
for achieving optimality for a convex optimization problem. 
We intend to use a primal-dual interior point algorithm for 
problem (1221 . which uses Mehrothra's predictor-corrector 
scheme [21]. The KKT optimality conditions which result 
from (l22l are: 



Hz 



C T v + G T \ = 
Cz - c = 
Gz - b + s = 
AS = Q 
A > 0,s > 0, 



where s £ R q are slack variables, v e M" M and A € K q 
are the Lagrange multipliers and S — diag(s), A = diag(A) 
are diagonal matrices formed from the slack variables and 
respective multipliers. These conditions lead to the following 
Newton system (see [1] for more details): 



(25) 



We can eliminate As by using As = —A 1 (r s + SAX). 



H 


C T 


G T 





C 











G 








/ 








s 


A 





~Az~ 




r z 




Av 




r v 




AX 




r\ 




As 




r s _ 



Furthermore, by reducing AA = S 1 A (r\ 
we obtain the following system: 



where 
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(26) 



$ = H + G T S~ 1 AG 
r d = r z + G T S- 1 Ar x -G T S- 1 r s . 

Next, we form the Schur complement of the matrix in 
so as to obtain the final system of equations: 



(27) 
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Solving ( I28t would normally employ the computation of 
Y, which may appear to be overwhelming given the large 
dimensions of Y and the fact that it requires an inversion of 
$. However, due to the way in which matrix Y is formed in 
d29t . we show that we can compute its Cholesky factorization 
in an efficient and distributed manner, similar to the one 
found in [3] for one linear system. The matrix $ £ M nxn 
has a block-diagonal structure $ = diag($ l ), where the 
blocks $' £ K""'^" 1 are also block diagonal, with their 
first block of size m, x mi, the following N — 1 blocks of 
size rii x rii and the final block of size n,; x rij. Now, it 
can be observed that resulting matrix Y has the following 
block-tridiagonal structure: 
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where the matrix blocks are: 
Y n = C 11 ^ 1 )- 1 ^ 11 ) 2, 

yii _ ^i-l^i-lyl^i-l^T (32) 

+ (7"($ l )" 1 (C") T , Vi = 2...M 
Y i,i-i = c i > i - 1 {$ i - 1 )- 1 (C i - 1 > i - 1 ) T ,Vi = 2...M. 



First, we show how to compute efficiently in a distributed 
fashion the matrix Y. Note that inverting the block com- 
ponents of $ and then forming the block components of Y 
would be very inefficient. However, if we form the Cholesky 
factorization of <I> 1 = L l (L nT 



W 



we get: 

V" = c u (L i y T 

1 ^ji.i— 1 (x i—1 \ 



(l )~ 



(33) 
(34) 



where L 8 g ]jjA r n i xA'n i a j so block diagonal, so that the 
block components of Y are: 



r 11 = 



= w 



1 fxri— i.i— 1 \T 



T, V 4 >2 



The most efficient computation of V t ' l ~ 1 can be done by 
solving the following systems of matrix equations, where the 
matrices Lp with j = 0, . . . , N, are the diagonal elements 

m 



of L 1 



i ^ = (A l ) T , 
r = (Ef,Vj = l 

T 



.N-l 
.N — 1 



(35) 
(36) 
(37) 
(38) 



Equations d35ll-(l3~8ll can be efficiently solved by matrix 
forward substitution, due to the lower triangular form of Lj. 
The resulting matrix will take the following form: 
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To obtain W 1 ' 1-1 we solve the following series of matrix 
equations, by matrix forward substitution, considering the 
dense nature of B 1 ^ 1 and A M_1 : 



3+1,3 + 



1 ) T = (A*'«- 1 ) T .Vi = l. 



(39) 

,N-1, (40) 



where L l 1 are the diagonal elements of 1/ 1 . The resulting 
W' l ' l ~ matrices will have a block-diagonal structure. 

Second, the resulting structure of the Cholesky factoriza- 
tion of Y — LL T is as follows: 
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L" 1 can be obtained from the following: 

-lXf T ll\T _ vll 

V' % -\L 
L U (L H ) T 



L ll {L ll Y =Y' 

-i.i-lt T i-l,i-l\T 



Y l 



(41) 

= Y 1 '^ 1 (42) 
L i ' i - 1 (L i ' i - 1 ) T , Vi > 2. (43) 



Note that the matrices Y %% have a block tridiagonal struc- 
ture, but L % ' % ~ x are usually dense so that there is no special 
structure in the terms Y n — Z/' I_1 (L M_1 ) T '. Therefore, the 
Cholesky factorization of these matrices is computationally 
demanding. 

V. Discussion on implementation 

The most important aspect of the algorithm previously 
presented is that it can be implemented in a distributed 
manner, between the M subsystems. The Cholesky fac- 
torization of Y clearly dominates the system of equations 
(l28l i when it comes to computing cost. By computing the 
matrices V M and W l,% ~ x , the inversion of $ can be avoided, 
and they can further be used in (l30l > and (I3H to calculate 
the respective residuals. The factorization of Y is also the 
most complex when it comes to the communication between 
subsystems, requiring the back and forth transmission of 
matrices between subsystems. Also, due to the structure of Y , 
the factorization cannot be done in parallel and is achieved in 
a sequential manner. For subsystem i, with i = 2, • • • , M — 1 
the following steps are required for obtaining L il and L 1 ' 1 ^ 1 : 

1) Compute $j = LjLf 

2) Send 1/ to subsystem i + 1 

3) Receive L l_1 from subsystem i — 1 

4) Compute V u : solve (fJSJ to © 

5) Send V %% to subsystem i + 1 

6) Compute W 4 '* -1 : solve ((39) and gOj 

7) Compute Y n , receive V l ~ l ' % ~ 1 from subsystem i — 1 

8) Compute Y % ' % ~ 1 , receive L t-1,t-1 f rom subsystem i — 1 

9) Compute L 1 ^ 1 from (@3 

10) Compute L l% from d43l . send L 1 ' 1 to subsystem i 

The number of flops for computing the Cholesky factoriza- 
tion of Y by each subsystem are provided in Table [0 

TABLE I 

Number of flops for computing local components of 
the Cholesky factorization of Y 



Operation 



\fi = 2 . . . M. The block components L l% and 



Factor: = I/(I/) T 
Solve: (35) 
Solve: (36), (37) 

Solve (38) 
Solve (39} 
Solve (40) 
Compute: K" 
Compute: yM-i 

Compute: L" 
Compute: L !,l_1 



Number ot Hops (approximate) 



3 1 3 

mm? 

2(N - l)n 4n 2 

nf 
3 

(TV - lKnf.j 
ATn?(n» + nj+ni_i +2) 
N(niiii-iai-i) 

JV 3 nJ 
3 



It can be observed that the matrices transmitted back and 
forth are very sparse, with a known block structure such 
that the only data required to be transmitted are these 
comprising blocks. Also, these blocks are transmitted only to 
neighboring subsystems, such that the transmission of data 
is localized. 

Note that the cost of computing matrices L' a and 
is cubic in N but linear in M overall, given the choice 
of z. Also, computations can be done sequentially and 
exchange of information is only between neighbors. If we 
would rearrange z by the prediction horizon, instead of by 
subsystems, then the dominating cost for computing these 
matrices would be linear in N overall and cubic in M locally, 

i.e of order ^ — >= * for L 11 . However, this would imply 
that every subsystem has knowledge of the dynamics of all 
other subsystems, and as a result computations would require 
all-to-all transmission of data between subsystems. Thus, the 
efficient choice of z given a physical leader-follower system 
involves the imposed prediction horizon N, the number 
of the subsystems M and possible transmission limitations 
between subsystems. 

VI. Conclusions 

In this paper we have showed that by restructuring certain 
MPC problems for large-scale systems we can reduce the 
computational cost of implementing an interior point algo- 
rithm foe solving such problems. An analysis for obtaining 
a stabilizing linear control law from a distributed viewpoint 
has been made. By combining several recent results, we have 
proved that the online computation of MPC control laws for 
some special classes of large scale systems can be carried 
out with increased speed through a reduction of the number 
of required flops. This, in combination with ever-increasing 
distributed computing power that can be used for distributed 
computation of an MPC law suggests us that MPC can be 
used now in many large-scale applications where it has not 
been considered applicable before. 

Further details regarding the efficient transmission of data 
between subsystems and the implementation results for the 
interior point method presented are omitted for lack of space. 
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